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Abstract 

The return distributions of the coherent noise model are studied for the system 
size independent case. It is shown that, in this case, these distributions are in the 
shape of g-Gaussians, which are the standard distributions obtained in nonexten- 
sive statistical mechanics. Moreover, an exact relation connecting the exponent r 
of avalanche size distribution and the q value of appropriate g-Gaussian has been 
obtained as g = (r + 2)/r. Making use of this relation one can easily deter- 
mine the q parameter values of the appropriate g-Gaussians a priori from one of 
the well-known exponents of the system. Since the coherent noise model has the 
advantage of producing different r values by varying a model parameter a, clear 
numerical evidences on the validity of the proposed relation have been achieved for 
different cases. Finally, the effect of the system size has also been analyzed and an 
analytical expression has been proposed, which is corroborated by the numerical 
results. 
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1 Introduction 



Throughout the last two decades the interest in extended dynamical systems has ex- 
perienced a steady increase. These systems exhibit avalanches of activity whose size 
distributions are of power-law type. Although there is not a unique nor unified theory 
which totally explains all the features of these complex systems, there exist several 
known mechanisms producing power-law behavior. One of the most popular and well- 
studied mechanisms is that of self-organized criticality (SOC) introduced by Bak, Tang 
and Wiesenfeld Many physical systems and models have shown to exhibit SOC 
llJ]. The most important feature of all these systems is that the entire system is under 
the influence of a small local driving force, which makes the system evolve towards a 
critical stationary state having no characteristic spatiotemporal scale, without invoking 
a fine-tuning of any parameter On the other hand, SOC is not the only mechanism 
causing power-law correlations that appear in a nonequilibrium steady state. Another 
simple and robust mechanism exhibiting the same feature in the absence of criticality 
is the coherent noise model (CNM) miU. The CNM is based on the notion of an ex- 
ternal stress acting coherently onto all agents of the system without having any direct 
interaction with agents. Therefore, the model does not exhibit criticality, but it still 
gives a power-law distribution of event sizes (avalanches). 

Recently, it was presented an analysis method to interpret SOC behavior in the lim- 
ited number of earthquakes from the World and California catalogs by making use of 
the return distributions (i.e., distributions of the avalanche size differences at subse- 
quent time steps) [5\. In their work Caruso et al obtained the first evidence that the 
return distributions seem to have the form of q-Gaussians, standard distributions ap- 
pearing naturally in the context of nonextensive statistical mechanics l6j|7j. Based on 
the assumption that there is no correlation between the size of two events, they were 
also able to propose a relation between the exponent r of the avalanche size distribution 
and the q value of the appropriate g-Gaussian as 



which is rather important since it makes the q parameter determined a priori and there- 
fore it acquits q of becoming a fitting parameter. The only little drawback of their 
work was that the number of data taken from the catalogs is not sufficiently large to 
obtain a very precise r exponent and also clear return distributions with well-defined 
tails (which is important in order to verify how good the distribution approaches a q- 
Gaussian). Consequently, Eq. ([T]) could not be rigorously tested until a very recent 
effort by Bakar and Tirnakli in |8|, where the same analysis was made using a sim- 
ple SOC model known as the Ehrenfest dog-flea model in the literature [9J (see also 
ifTOl fTTI ). Thanks to the simplicity of the dog-flea model, it was possible to achieve 
extensive simulations with very large system sizes (up to lO*^) and also very large num- 
ber of data elements (up to 2 x 10^). Accordingly, from these extensive simulations, it 
was obtained a value of t — 1.517, which is in accordance with the "mean-field" ex- 
ponent 3/2 determined in several problems lfT2l [131 [T4l . Thence the q value of return 
distributions was deduced a priori from Eq. ([T]l- 

In this work, we plod along this way by setting forth the following points: (i) first, 
we will obtain an exact relation between r exponent of the avalanche size distribution 
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and the q value of the appropriate g-Gaussian without resorting to any assumption 
and compare it to Caruso et al relation given in Eq. ([T]), (ii) since the CNM has the 
advantage of producing different r values by varying a model parameter crQ we now 
have the opportunity to test the validity of our exact relation (and also the Caruso et 
al relation) not only for one case but for various cases, (iii) since the corresponding 
return distributions are expected to converge to the g-Gaussian as the system size goes 
to infinity, the effect of finite system size is also important and we shall try to analyze 
this effect proposing an analytical expression, (iv) and finally since this model is not a 
SOC model, our results also give us the possibility of checking the generality of this 
behavior observed so far in SOC models. 

2 The coherent noise model 

Let us start by introducing the CNM. It is a system of N agents, each one having a 
threshold Xi against an external stress 77. The threshold levels and the external stress are 
randomly chosen from probability distributions pthresh (x) and Pstress (?7), respectively. 
Throughout our simulations we use the exponential distribution for the external stress, 
namely, Pstress = (l/o") exp(— t^/ct) and the uniform distribution (0 < x < 1) for 
Pthresh(x). The dynamics of the model is very simple: (i) generate a random stress 77 
from Pstress (v) replace all agents with Xi < rjhy new agents with new threshold 
drawn from pthresh{x), (ii) choose a small fraction f of N agents and assign them 
new thresholds drawn again from pthresh {x), (iii) repeat the first step for the next time 
step. The model can be described in the form of a two step-master equation that we 
present in the appendix. The number of agents replaced in the first step of the dynamics 
determines the event size s for this model. Although the CNM has been introduced for 
analyzing biological extinctions ||3], it has then been adopted as a very simple mean 
field model for earthquakes even though no geometric configuration space is introduced 
in the model [4 |. It is shown that the model obeys the Omori law for the temporal decay 
pattern of aftershocks lITSl . exhibits aging phenomena lfT6l and power-law sensitivity 
to initial conditions f\n\. 

3 Avalanche Size and Return Distributions 
3.1 Size independent case 

As pointed out in |4 |, there is advantage in choosing the uniform distribution (0 < x < 
1) for the thresholds of the CNM agents seeing that the model can be simulated in the 
N 00 limit using a fast algorithm which acts directly on the threshold distribution 
instead of acting on the agents of the system. This enables us to obtain the avalanche 
size distribution of the model as being independent of the system size. The distri- 
bution P{s) is expected to be a power-law over many decades until the s values reach 
a particular point s ^ cr, thereafter it falls off exponentially. From our point of view 

'in the dog-flea model there is only one available value of t since the only parameter is the number of 
fleas. 
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this is rather important since it means that if we measure the avalanche size exponent 
T using the region s < cr, then we must use this t value to predict a priori the q value 
of the q-Gaussian that the return distribution is expected to converge in the entire re- 
gion without any deterioration (not only in the central part but also in the tails). The 
results obtained for the avalanche size distributions of three representative cases with 
a — 0.01, a — 0.05 and a = 0.065 are given in the left column of Fig. 1. 

Since each case with different a values has a different size exponent t, this al- 
lows us to check the validity of Caruso et al relation given in Eq. ([T) or any other 
equation relating r values to the q values of the appropriate g-Gaussians. From the 
master-equation of the CNM is theoretically possible to compute the probability of s 
and bringing to bear standard techniques 1 13| to obtain the return distribution. How- 
ever, its level of complexity turns out the solution almost analytically impossible or its 
(asymptotic) behavior deeply unclear as it happens in several other problems of this 
class 1 18 1. Regardless, we are in the position where we can propose an exact relation 
for the return distribution P (As) bringing into play no other assumption than the dis- 
tribution of avalanche sizes, where As is the difference between two consecutive event 
sizes, i.e.. As = s{t + 1) — s{t). Let us mathematically define the avalanche size 
distribution, 

p(s) oc (e + s)-", (r>l), (2) 

with e being a constant value describing the asymptotic limit s — ?> 0. The process of 
avalanches is completely Markovian (independent) and therefore the probability of the 
difference of sizes As is 

/•oo /"OO 

P(As) = / / p{s)p{s')S{As-{s~s')) ds'ds. 
Jo Jo 

(e + sy^ (e + As + s)"^ 9 (As + s) ds, 

where 0(. . .) is the Heaviside step function and s' denotes the previous avalanche size. 
Making use of |19| and attending to the symmetric nature of P (As) we can explicit 
the negative branch, 

F(As) = |As|i-2^(B[^,l-r,l-r] (-1)^+C(t)), (3) 

where C (r) is a coefficient only depending on t and related to the convolution of 
very large values of s with very large values of — s' yielding a |As|^ '^^ dependence 
0. Thus, the distribution is mainly described by the product of the isolated factor by 
the incomplete Beta function . .]. Applying the asymptotic behavior a: — of 
B[x,a,b] Ii20il we finally get, 

P(As)^|Asr", (As»l). 

Taking note of the g-Gaussian distribution, 

P(As) = P(0) [1 + ^{q - l)As2] '^^'^"^ , (4) 



This can be flatly checked out performing the calculation with e = 0. 
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we straightforwardly obtain 



This relation is slightly different from the approximate relation presented in 15] as can 
be seen in Fig. 2. For t ^ oo, both relations approach q = 1 and they are almost 
identical except in the region where r values are smaller than 1.5. Moreover, only the 
relation (|5]l correctly achieves q = 3 value when t = 1. These values define the limits 
of the domain of each parameter so that the distributions dU and are normalizable. 
The approximate relation Eq. ([T]i does not fulfill this condition as g (r = 1) > 3. It 
should be noted that, since this discrepancy is only meaningful for r < 1.5, the ap- 
proximate relation predicts the q values with a |0.01| difference from the exact one, 
which are also acceptable for all the cases we present. 

We are now ready to proceed analyzing the return distributions. The centered re- 
turns are given in terms of variable x 

a; = As -(As) , (6) 

where (...) represents the mean value of a given data set. As can be seen from the 
right column of Fig. 1, in our simulations we generated the return distributions of the 
three representative cases of the CNM in order to check the validity of the relation (|5]l. 
In each case, an extremely large number of events (2 x 10^) has been used to build 
the numerical distribution, namely, the central part and tails. It is clear that the return 
distribution (green dots) can by no means be approached by a Gaussian. They actually 
exhibit fat tails which agree with (/-Gaussians Eq. (|4|i where /3 characterizes the width 
of the distribution and q is the parameter which should be determined directly from 
Eq. Q a priori and therefore is no longer a fitting parameter. In each panel on the 
right column of Fig. 1, the dashed black lines represent the appropriate (/-Gaussian 
with the q value obtained from Eq. (|5]). Perfect agreement with the data can be easily 
appreciated not only for the tails but also for the intermediate and the very central part 
as it is demonstrated in the insets. 



3.2 Size dependent case 

Although we might think that the size independent (i.e., infinite size) case would be 
enough for such an analysis, we believe that it is still instructive to look also at the size 
dependent case at least from two different perspectives: (i) we can check how the size 
of the system affects the shape of the return distributions and whether the tendency is 
consistent with the infinite size case as the size of the system increases, (ii) unlike the 
CNM, generic size independent cases cannot be achieved for such systems and thus the 
only possibility is to always analyze the size dependent case. 

It is very easy to implement the size dependent algorithm for the CNM. We just 
need to apply the previously described steps of the dynamics to a system of N agents. 
As N increases, this algorithm clearly slows down and for the same number of events 
(10^) the larger value of N that we can simulate in a reasonable time is 20000. In 
Fig. 3(a) the behavior of the avalanche size distribution is given for a — 0.05 case for 
various values. It is clearly seen that the power-law regime is always followed by 
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an exponential decay of all the curves and this decay is postponed to larger sizes as 
N increases. For each N case, we estimate the r value using the standard regression 
method in the region before the exponential decay (we determine the size of this inter- 
val by looking at the regression coefficient to become always more than 0.9997 in each 
case). Therefore, we should expect that the exponential decay part would tamper with 
the g-Gaussian behavior of the return distributions and this meddling must diminish as 
N gets larger and larger, which is in fact observed in Fig. 3(b) for the return distribu- 
tions of four representative N values. When N values are very small, avalanche size 
distribution has a very short power-law region and the exponential decay part dom- 
inates, which simply causes the return distributions to deviate immediately from the 
g-Gaussian shape. As N increases, return distributions start approaching the thermo- 
dynamic limit (dotted black line), which is a full g-Gaussian with q = 2.09, yielding 
better and better from the central part to the tails, i.e., as the expected scale- free regime 
sets in. 

In order to explain this gradual approach to q-Gaussians when finite-size effects are 
present, let us try to develop a simple mathematical model by considering the differen- 
tial equation 



dy 
d(a;2) 



-aru'' - {aq - ar)y'' (aq > ar > 0; q > r; y(0) = 1) 



(7) 



This equation has very interesting and different solutions depending on the choice of r 
and q values (see refs. 1171 12111221 ). but for our purpose, let us concentrate on case r — 1 
and q > I, whose solution is given by 



y = 



1/(1-9) 



(8) 



If oi = 0, then the solution coincides with the q-Gaussian, whereas if Uq = ai 
(which means that q — 1), the solution turns out to be the Gaussian. On the other 
hand, between these two extremes, namely if a, > ai > and g > 1, we obtain a 



crossover between them 
g-Gaussian, y ^ [l — (1 



Specifically, for {q 

,211/(1-9) 



1) fli x"^ ^ 1, Eq. ^ approaches a 
q) aq a;^] ^^'^ Our results, which are depicted in Fig. 3, 
show that the small values of ai imply that the (/-Gaussian form is valid up to rather 
large values of x. On the other hand, for {q — l)ai x^ 3> 1, the exponential outnumbers 
the remaining terms leading to the Gaussian behaviour. 



1/(1-9) 



exp ( — ai X 



The approximate dependence of Eq. (O can thus be split into different regions defined 
by three values of x. Namely the first value is 



W 
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(W [. . .] is the Lambert W function ll20ll ) whence the curve assumes a power-law de- 
pendence described by the exponent 2/{q — 1) that persists up to 



when the it starts being perturbed by the Gaussian dependence. Last, there is the final 
convergence to the Gaussian functional form which occurs at 



This crossover seems to coincide with the behavior of the return distributions of 
the N dependent cases as plotted with dashed black lines on top of each curve in 
Fig. 3(b). This behavior simply reveals that the longer the power-law regime per- 
sists for avalanche size distribution, the better the appropriate q-Gaussian dominates 
in the return distribution. Finally, as — > oo, the power-law regime prevails for the 
avalanche size distribution giving forth a return distribution following the appropriate 
g-Gaussian for the entire region. 



In this work, we have studied the behavior of the return distributions for the CNM by 
directly simulating the size independent case. By means of extensive simulations, it is 
clearly shown that these distributions converge to g-Gaussians with appropriate q val- 
ues which are deduced a priori from the exact relation Q that we developed here. It is 
worth noting that although the q-Gaussian description is actually an analytical approx- 
imation the result provides for an understandable depiction of the distribution, which 
hardly occurs when we keep a special functions representation, with no fundamental 
accuracy lost. This relation makes the q parameter be related to one of the well-known 
exponents (avalanche size exponent r) of such complex systems and therefore it res- 
cues q from being a fitting parameter in this analysis. Moreover, since the model pa- 
rameter a allows us to obtain different r values, we were able to check this behavior for 
various cases. These results clearly imply that the observed behavior is not restricted 
to self-organized critical models, but instead it seems to be a rather generic feature pre- 
sented by many complex systems which exhibit asymptotic power-law distribution of 
avalanche sizes. 

We have also investigated the finite-size effect by simulating directly the model 
dynamics and found that the convergence to appropriate g-Gaussian starts from the 
central part and gradually evolves towards the tails as the system size increases. This 
is in complete agreement with the gradual extension of the power-law regime in the 
avalanche size distribution before the appearance of the exponential decay due to finite- 
size of the system. These results corroborate the analysis of size independent case since 
it is clearly seen that, as N oo, curves of return distributions for size dependent case 
converge to the one comes from the size independent case. 





4 Conclusion 
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Finally it should be noted that, since it is generically extremely difficult (if not 
impossible) to achieve the size independent case for such complex systems, the size 
dependent case has its particular importance. Therefore, although the return distribu- 
tions appear to be g-Gaussians for the entire region in the thermodynamic limit, we 
have tried to propose a mathematical model in order to explain the behavior of return 
distributions for the size dependent case. 
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A The CNM master equation 

The dynamics of the CNM can be described according to the probability of having 
n agents in the system that at time t present a critical value up to x, P„ {x,t). In 
conformity with step 1 we can write the master-equation, 

n 

Pn{x,t) = P„(a:,i-l)W'„^„+^P„_;(a;,t-l)M^„_j^„ 

1=1 

N — n n 
+ J2 Pn+l t - 1) Wn+l^n - ^ P„ (x, t - 1) Wn^n^l 
1=1 1=1 
N-n 

- ^ P„(a:,t-l)W„^„+, (9) 
1=1 

with the probability transitions W given by 

Wr^^n = F' (x) + V / x"' - pthresh (v) 6 (m - M,,^) dyy, (10) 

where the first term on the rhs comes from the case i] > x and the second one otherwise. 
The inverse cumulative probability F' (x) = 1 — F (x) — 1 — Jq Pthresh (z) rfz^and 

Mx.rj = '^iLi '^[v ~ Xi {t)] Q[x — ij] means the number of agents with critical value 
below X and 77. The following elements are 

Wn^l^n = Wn^n+l = 0, (11) 



'For our case, i.e., Pthreshi^) = cr ^ exp[— z/ct] implies F(z) = 1 — exp[— z/ct]. 
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Wn+l^n - ( " + M (1 - X)' F' {X) + 



n 



ri+l rx / \ 



- 1 



a;"-' (1 - xy F' (x) + 

This corresponds to a matrix with vanishing elements below the diagonal. From 
these relations is then possible to spell out the occurrence of an avalanche of size s 

Z^«=o Jo dx ^^n^7i ax <^ S - u 



Pis) = { ^N-s rl dP„+4x,t-i) , (14) 

= s ^ 



2^n=0 Jo dx ''yn+s^n ax-^- 

rl dP„(x.t-l) , 



which is numerically well described by the power-law (|2|l with a small value of e. 
Regarding step 2 the master equation is abstractly pretty much the same, 

n N-n 
Pn {X, t+l)^ Pn {X, t) Wn^n + E Pn-l (x, t) Wn^l^n + E Pn+l (x, t) W^+l^u 

1 = 1 1 = 1 

n N—n 
~Y.Pn {X, t) W^^n^l ~ Pn {x, t) Wn^n+l , (15) 

1=1 1=1 

with the probability transition matrix is given by 

,V " \( N~n \f fN \( nyf^ 







n 



{l-p)fN 

{l-x)[l-j^)\ Q[n-pfN]Q[N-n-{\-p)fN]dp (16) 

where p is used to define the subfraction of agents, pJN, whose critical value before 
updating was less than x. 



pfN ) i {l-p)fN ) i pfN + l 



n 



■PfN+l n _ ^^(l-p)/JV-^ 



N ) V N 

Q[n-l- pfN] e[N - (n-l) - (1- p) fN] dp (17) 
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M/ _ fW n + l \ f N - (n + l) \ f fN \ 

Wn+l^n- ]^ U/^ J I (l-p)/iV ) \ pfN-l ) ^ 

/ , i \ pf N / , ,\ n-p)fN 

e [n + / - pfN] e[N-{n + l)-{l-p) fN] dp (18) 



n— >n— i 



ViV 





PfN 



n \ f N-n \ f fN 
PfN ) \ {l-p)fN ) \ pfN- 

/ n\(i-p)/^ 
i}-N) 

Q[n-pfN] Q[N -n-{l- p)fN]dp; 



(19) 



n— )-nH-/ — 



nW AT-n \/ fN 
PfN ) \ {l-p)fN ) \ pfN + l 



(n) v-n) ^"'"^'a-^y 

e[n- pfN] e [iV - n - (1 - p) fN] dp. (20) 
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Figure 1 : Left column: Avalanche size distributions of three representative values of a. 
For each case, the r value is calculated using standard regression method for the region 
s < (J. Right column: Return distributions for the same three cases. Two zooms of the 
central part are given in the insets for better visualization. For each case, / = 10^^ 
and 2 x 10^ experiments are generated. 
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exact relation in Eq.(5) 

Caruso et al. relation in Eq.(l) 




Figure 2: Comparison of our exact relation (Eq. (|5]l) with that of Caruso et al (Eq. ([T}). 
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Figure 3: (a) Avalanche size distributions for N dependent case, (b) Corresponding 
return distributions of the same N values. 



14 



